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Abstract. - A "sandpile" cellular automaton achieves complex temporal correlations, like 
a 1// spectrum, if the position where it is perturbed diffuses slowly rather than changing 
completely at random, showing that the spatial correlations of the driving are deeply related 
to the intermittent activity. Hence, recent arguments excluding the relevance of self-organized 
criticality in seismicity and in other contexts are inaccurate. As a toy model of single fault 
evolution, and despite of its simplicity, our automaton uniquely reproduces the scaling form of 
the broad distributions of waiting times between earthquakes. 



Many complex natural phenomena exhibit quiet periods with a slow loading of the sys- 
tem separated by events in which a rapid evolution takes place [1,2]. Some examples are 
earthquakes [3,4], solar flares [5], creep experiments with cellular glasses [6] and even with 
piles of rice [7], and rain falls [8]. In these and in other natural and social systems one finds 
that the intensity of the events and the spatial and temporal scales separating them may 
vary over broad ranges and often are scale-invariant. The similar phenomenologies and the 
nonequilibrium nature suggest that a common mechanism is emerging in different systems. 
An interesting candidate is self-organized criticality (SOC) [1,2,9], a theory essentially fo- 
cussing on the occurrence of fast avalanches of nonlinearities, which have been well visualized 
in simple cellular automata (CA) called "sandpiles" . 

CA and their probabilistic versions are model-systems that have been used quite exten- 
sively in the field of complex systems [10,11]. Good reasons include the possibility of reli- 
able simulations and the conceptual simplicity of the dynamical rules. Furthermore, one has 
strong indications that the microscopic details of nature are not all-important in deciding cer- 
tain macroscopic features. Symmetry properties, conservation laws and considerations on the 
right scale of description matter much more. CA incorporate these essential ingredients and 
express in the most strongest form we know today how emergent behavior can be very rich, 
varied and complex even if based on few simple dynamical rules on a discrete architecture. 

A feature distinguishing several sandpile models from other CA is the scale- free distribution 
of their avalanche sizes. This resembles the distribution of energy released by earthquakes, 
for example. On the other hand, temporal correlations between avalanches have been far 
less investigated, until Boffetta et al. noticed that realistic correlations are not found in time 
scries of some sandpiles [12]. Indeed, their avalanches have interoccurrence times distributed 
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exponentially, and the time series appears as a Poisson process of uncorrelated events. This 
fact was used to argue that SOC cannot be an interpretation of solar flares activity [12] and, 
more recently, that SOC is not related to earthquakes [13]. 

After the critique to SOC by Boffetta et al. [12] there has been a wave of investigations on 
the temporal properties of sandpiles and other SOC models. Nowadays we know that there 
are models with interesting temporal correlations, like the 1// decay in power spectra [14-16], 
or like power-law tails in distributions of waiting times between events [17-22], foreshocks and 
aftershocks like for earthquakes [23] or even features typical of turbulence [24,25]. 

Correlated drivings can give rise to non-exponential waiting time distributions between 
events [18], which e.g. appeared in models of solar activity. In one case the strength of the 
driving perturbation in the Lu-Hamilton model was slowly varied at random [17]. In another 
model of "Ellerman bombs" [19], directed percolation was one mechanism increasing the 
number of driving points. These examples further confirm and illustrate the coexistence of 
criticality and complex time correlations but it is likely that the correlations intrinsic in these 
driving mechanisms are simply inherited and show up directly in the output of the automaton. 

In this Letter we stress that a completely random driving of sandpiles is not expected to 
be a realistic feature but rather an artificial feature put by hand. In other words, random 
drivings can be an arbitrary, a priori limitation of the automata. Usually sandpiles are driven 
either at completely random positions or from a fixed single site, but none of these cases 
reflects the forms of slow loading in real systems, like the crust of the Earth or the solar 
corona. In seismicity, for example, one observes a rich scaling picture for the waiting time 
distributions [3,4], involving location, temporal occurrence, and magnitude of events. In order 
to have satisfactory CA models of earthquakes, it is therefore desirable to remove excessive 
randomness from the properties of their driving. 

Below, we show how extremely simple spatial correlations in the input drive can generate a 
complex time series as an output. Thus, the focus is on the spatial properties of the nucleating 
point of the avalanches. As far as we know, we provide the first example of 1// noise in 
the avalanche activity of a non-running, classical sandpile. Furthermore, we observe broad 
distributions of waiting times between events also of small sizes. Most importantly, these 
distributions display a scaling behavior with intensity thresholds, analogous to those recently 
found for earthquakes [4] and for solar flares [5] . This feature makes at present this new model 
unique, and corroborates its interpretation in geophysical terms. 

The dynamics of the model strongly resembles the one of the sandpile originally discussed 
by Kadanoff et al. [26]: each site i of a linear lattice of side L holds an integer number hi. 
In the sandpile jargon, hi is normally called "number of grains" , but this distribution as well 
may be thought as the potential energy profile in an object mainly extended in one dimension, 
like a fault. Stable configurations fulfill a local stability condition on the discrete gradient 
between every pair (i, j) of nearest neighbors, 

hi -hj <H (1) 

where constant H is a threshold. A local instability not fulfilling Q is resolved by a toppling, 
which consists in moving a grains from i to j, i.e., hi — > hi — a and hj — > hj + a. In the 
geophysical interpretation, the toppling thus models a redistribution of energy that takes place 
after some local elastic threshold within a fault is overtaken. 

Stress increase appears as a local increase of potential energy during time step t at position 
i', causing he (t) = — 1) + 1. This eventually leads to a violation of JD between site i' and 
one or more of its nearest neighbors j. Thus, a toppling takes place between each unstable 
pair and j in turn may also become unstable. One then iterates this procedure by 
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Fig. 1 - (a) Snapshot of a sandpile profile (L — 256) around its average height. The arrow points to 
the position i' where grain addition was taking place, (b) Example of time series of the avalanche 
sizes for a sandpile with L = 2048. 



making all topplings in parallel as long as some pairs of sites violating (JJJ are present [27]. 
The whole set of updates giving rise to a new stable configuration represents an avalanche 
(earthquake) at time t, with an area a equal to the number of sites that toppled at least 
once and a size s (energy released) equal to the total number of toppling. As in the standard 
sandpile scenario, the long time scales of the driving are thus completely separated from the 
ones of the avalanches. This is the appropriate limit for models of earthquakes, while mixing 
of time scales is present in solar activity [5], and other models of SOC are possible [22]. 

A novel feature distinguishing the model proposed here from previous ones is that the 
site where a grain is added coincides with the position of a random walk: i'(t) with equal 
probability is drawn from one of the nearest neighbors of i'{t — 1). This feature again has a 
counterpart in seismicity, where it is known that epicenters are spatially clustered and that 
aftershocks slowly diffuse after large events [28]. 

The particular driving introduced here makes the model critical also with periodic bound- 
ary conditions (BC) [29]. We adopted periodic BC because there cannot be effects related to 
a fluctuating distance between a predefined boundary and the point where grains are added. 
For our purposes it is enough to show results on the sandpile in one dimension, with H = 4 
and a = 2. According to our data, the versions with different H or a give similar results. The 
same remark applies to choices of different mobility; as long as the walker diffuses slowly, no 
change occurs in the basic aspects of what follows. 

A typical sandpile profile is plotted in Fig. da) ■ Configurations like that are reached after 
an appropriate period of fueling of an empty lattice. That period was then not considered in 
the following statistical analysis. The sandpile profiles observed with periodic BC usually are 
an alternated sequence of uphills and downhills, thus alternating local minima and maxima. 
Each patch with the same trend is like a profile of a classical sandpile with open BC, which 
oscillates around a typical slope. However, grain addition and avalanches make patches to 
dynamically evolve, merging them or splitting them into shorter pieces. 

A time series of the avalanche sizes in the stationary regime is shown in Fig. QJb) . It has 
an intermittent behavior with bursts of activity and temporal clustering of large avalanches, 
features due to time correlations, as shown below. The probibility density of sizes for several 
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Fig. 2 - Probibility density of the area of avalanches, for lattices with L — 256, 512, . . ., 4096, from 
left to right. The small secondary tail of P(a) is coming from rare avalanches with a > L/2 involving 
grains falling onto two sides of a mountain. Insets: (a) Size distributions for the same set of L (curves 
are smoothed), (b) Collapse of Pl(o), according to Eq. with t = D = 1, onto a scaling function 
F. 



L's, -Pl(s), are shown in Inset (a) of Fig. [21 and display multiscaling [30] while the area 
probability density Pl (a) (Fig. |2J) obey simple finite size scaling 



as confirmed by a collapse onto a single scaling function F(x) in the plot of aPi(a) vs a/L 
(see the Inset (b) of Fig. The model thus has a critical behavior, with r ~ D ~ 1 [31]. 

Complex temporal correlations and long memory are revealed by the nontrivial form of 
the power spectrum S(f) of the time series of avalanche sizes s(t) (with s(t) = if at time 
t there were no topplings). Fig. Ufa) shows that at quite low frequencies S(f) ~ l// 7 with 
7 = 1.00(1). In Fig. E^b) we plot S(f) vs fL to show that the 1// region becomes broader 
when L in increased: it is bounded on the right by a bump at / ~ l/£, while the crossover 
frequency to a flat spectrum at very low / scales approximately as 1/L 2 . Almost equal spectra 
are observed for the time series of avalanche areas, confirming that the basic origin of that 1 // 
spectrum is in the temporal correlations of events rather than in their actual intensities. We 
remind that a 1// spectrum is observed in many natural phenomena [32], and it is regarded 
as one of the most complex forms of time correlation. It is remarkable that the long-range 
memory in the model arises by the non-trivial, self-organized structure encoded in the sandpile 
profile, in agreement with a basic idea of SOC theory, namely that noncquilibrium extended 
systems can self-organize to complex regimes even if their dynamical rules are simple. 

The other form of temporal correlation concerns the waiting time distribution between 
events which has a power-law tail, rather than an exponential decay as for uncorrelated events. 
The waiting time distributions P s (t w ) between avalanches of size > s in a sandpile with 
L = 4096 are plotted in Fig.^Ja) for some thresholds s ranging from 2 8 = 64 to 2 17 = 131072. 
Already for a small fixed s, the distributions have a power-law tail, which does not disappear 
when larger L are considered. In the Inset of Fig.^Ja), for example, one can see that P s =64(t w ) 
has the same power-law tail for L = 256 and L = 4096. Interestingly, in fact when s > L/2 
(steepest part of Pl(s), see Inset (a) of Fig. |2J) in P s (t w ) there are two regimes with two 
different power law decays, as one observes in the analysis of waiting times of earthquakes [4] 
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Fig. 3 - Log-log plot of the power spectrum of the time series of avalanche sizes, for L = 64, 128, 
256, 512, and 1024, (a) vs the frequency and (b) vs fL. In (a) curves have been offset half decade 
from each other, and a straight line indicates a 1// scaling. For the estimates of S(f), we used time 
windows with 2 21 data. 
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Fig. 4 - (a) Log-log plot of distributions of waiting times between events with size > s, for L — 4096 
and various s. Inset: the cases (L = 256, s — 64; •) and (L = 4096, s = 64; o) are compared, (b) 
Log-log plot of the same distributions, but rescaled with the rate of events > s. 



and for solar flares [5]. As in these cases, we find it appropriate to rescale times by multiplying 
them by rates of events > s, denoted as R(s) [33]. 

The rescaled times T — t w R(s) appear to be the natural ones, because P s (t w ) / R(s) vs 
t w R(s) collapse quite well onto a single scaling curve G(T) [see Fig. Efb)]. Since G{T)dT = 
P s {t w )dt Wl we see that G(T) is the probability density of T, not depending anymore on 
thresholds. Hence, a unifying scheme can be applied to all thresholds also in our sandpilc 
model, showing that there is a basic self-similar mechanism taking place. For short rescaled 
times G(T) - T^ 1 with ft = 0.8(1), while G(T) - with f3 2 = 2.0(2) for T > 1. (for 

earthquakes Corral found ft = 0.9(1) and ft = 2.2(1) [4]). The crossover between the two 
power-law regimes is around T* m 0.7. Thus, for t w much smaller than the average occurrence 
time t(s) = R(s)~ 1 of events with size > s (i.e. T <C 1), the correlations between events are 
qualitatively different from the ones at t w ^> t(s). 
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A mechanism giving rise to strong temporal correlations between avalanches should be 
their spatial overlap, since the memory of past activity is stored in the height profile of the 
sandpile. By dropping grains at a point slowly diffusing, as we do, a fresh part of the profile is 
almost always found, probably enhancing the correlation with past activity. In some classical 
models driven at random points, scale-free waiting time distributions are observed between 
sufficiently large avalanches [18], which indeed have a consistent probability to overlap with 
previous ones. On the other hand, that argument does not really explain why the 1// part of 
the spectrum of our sandpile extends to very low frequencies. 

We have mainly discussed SOC in CA, but a different class of SOC models also supports 
the view that space and time correlations are related to each other. We are referring to models 
with "extremal dynamics" . In these models each unit carries a continuum variable and the unit 
closer to instability is always the one relaxing (toppling). The Olami-Feder-Christensen (OFC) 
model of earthquakes [34] is representative of extremal dynamics. Whether criticality in the 
OFC model is only apparent or persists in the limit L — > oo is an open issue. Nevertheless, the 
OFC model has complex patterns in time [20, 23] and space: indeed, by its nature, the OFC 
map self-organizes also the position of the epicenters (driving points), which form sequences 
deeply linked to past activity [35]. Similarly, in the Bak-Sneppcn model of evolution, distances 
between extremal sites follow a Levy flight, and l// 7 spectra are found with < 7 < 1 [36]. 
Recently, Lippiello et al. introduced a hybrid sandpile/extremal dynamics model of seismicity 
that also display waiting time distributions with power-law tails [21]. 

In summary, we have studied a sandpile cellular automaton exhibiting self-organized crit- 
icality (even without open boundaries). A critical stationary state with correlated avalanches 
and intermittent transport takes place because the model is perturbed at a position slowly 
diffusing in space, a process leading to clustered nucleation points and motivated by the ob- 
servation of natural phenomena like earthquakes, where epicenters are correlated. Indeed, the 
distributions of waiting times between avalanches larger than given thresholds have striking 
and unique similarities with the ones found in real seismicity. We can thus reproduce some of 
the complexity of natural critical phenomena, involving jointly energetic, spatial and tempo- 
ral aspects. Our results suggest that the logic that led to the present model is right, namely 
that spatial correlations should not be arbitrarily removed by imposing a random drive in 
models of critical phenomena like sandpiles, which are by themselves already very simplified 
objects. The positive comparison of our results with earthquake activity leads to propose that 
a rugged landscape as depicted in Fig. ^a) may be related to the stress distribution along a 
fault. The absence of a single slope profile (as observed in sandpiles with open boundaries) im- 
plies that system-wide events are not always possible, with potential interesting implications 
about predictability of large events. 

* * * 
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